Start RStudio and open a new script file
Install required libraries in RStudio - if you do not have them already
install.packages(
c("sp","rgdal","tmap","classInt","RColorBrewer",
"ggplot2","leaflet", "ggmap"), dependencies=TRUE
)
Intro to working with geospatial data in R
Mapping geospatial data
Practice
Who are you?
Why are you here?
Open one of the tutorial files in a web browser
Slides r-geospatial-workshop-pt1.html
Tutorial r-geospatial-workshop-pt1-tutorial.html
Raw Code scripts/r-geospatial-workshop-pt1.R
Make sure you can cut and paste into RStudio
are data about locations on or near the surface of the Earth.
represent location more specifically with coordinates
46.130479, -117.134167
Coordinates only make sense when associated with a CRS!
Geographic Coordinates: Latitude and Longitude
Define:
the shape of the Earth
the origin (0,0 point)
the relationship between the system and the real world
the units
Because of variations in these, there are many geographic CRSs!
The World Geodetic System of 1984 is the most widely used geographic coordinate reference system.
WGS84 is the default CRS for most GIS software
Almost all longitude and latitude data are assumed to be WGS84 unless otherwise specified
Historical data much trickier
You can
dynamically determine spatial metrics like area, length, distance and direction
spatial relationships like intersects, inside, contains, etc
and
Spatial data is a broader term than geographic data.
Methods for working with spatial data are valid for geospatial data
All spatial data are encoded with some type of coordinate reference system
Geospatial data require the CRS to be a model of locations on the Earth
Vector and Raster Data
Points, lines and Polygons
Regular grid of cells (or pixels)
software that can import, create, store, edit, visualize and analyze geospatial data
represented as geometric data objects referenced to the surface of the earth via CRSs
with methods to operate on those representations
We call software for working with geospatial data GIS
Geographic Information System
This term is commonly associated with desktop software applications.
Desktop GIS - ArcGIS, QGIS
Spatial Databases - PostgreSQL/PostGIS
Web-based GIS - ArcGIS Online, CARTO
Software geospatial data support - Tableau
Programming languages with geospatial data support
R, Python, JavascriptYou already use R
Reproducibility
Free & Open Source
Strong support for geospatial data and analysis
Cutting edge
Make sure you have the packages we are going to use installed and loaded
Make sure you have set your working directory to the location of the workshop files
install.packages(
c("sp","rgdal","tmap","classInt","RColorBrewer",
"ggplot2","leaflet", "ggmap"), dependencies=TRUE
)
There are many approaches to and packages for working with geospatial data in R.
One approach is to keep it simple and store geospatial data in a data frame.
This approach is most common when
the data are point data in CSV files and
you want to map rather than spatially transform or analyze the data
sf_properties_25ksample.csv
San Francisco Open Data Portal https://data.sfgov.org
This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.
We are using a subset of these data as a proxy for home values.
sfhomes <- read.csv('data/sf_properties_25ksample.csv',
stringsAsFactors = FALSE)
# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
"NumBedrooms")]
Make sure your working directory is set to the folder where you downloaded the workshop files!
sfhomes <- read.csv('data/sf_properties_25ksample.csv',
stringsAsFactors = FALSE)
# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
"NumBedrooms")]
## YearBuilt totvalue AreaSquareFeet Neighborhood NumBedrooms
## 1 1923 1031391 2250 West of Twin Peaks 0
## 2 1909 775300 5830 Presidio Heights 3
## 3 1915 1116772 1350 Inner Richmond 2
## 4 1947 860130 1162 Sunset/Parkside 0
## 5 1981 322749 1463 Outer Richmond 3
class(sfhomes) # what is the data object type?
dim(sfhomes) # how many rows and columns
str(sfhomes) # display the structure of the object
head(sfhomes) # take a look at the first 10 records
summary(sfhomes) # explore the range of values
summary(sfhomes$totvalue) # explore the range of values for one column
hist(sfhomes$totvalue) # histogram for the totvalue column
Use the R base plot function to create a simple map
plot(sfhomes$lon, sfhomes$lat) # using base plot function
Use the R base plot function to create a simple map
plot(sfhomes$lon, sfhomes$lat) # using base plot function
ggplot2ggplot2Most widely used plotting library in R
Not specifically for geospatial data
But can be used to make fabulous maps
Great choice if you already know ggplot2
ggplot2Load the library
library(ggplot2)
ggplot2Basic map with ggplot
library(ggplot2)
ggplot() + geom_point(data=sfhomes, aes(lon,lat))
ggplot2Basic map with ggplot
ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1)
ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1) + coord_map()
Allows you to associate a map projection with geographic coord data.
Map Projection: mathematial transformation from curved to flat surface
A Projected CRS applies a map projection to a Geographic CRS
All introduce distortion,
in shape, area, distance, direction, or combo
the larger the area the greater the distortion
No one map projection best for all purposes
Selection depends on location, extent and purpose
totvalueData driven symbology
ggplot() + geom_point(data=sfhomes, aes(lon,lat, col=totvalue)) +
coord_map()
totvalueWhat’s happening here?
sfhomes_low2high <- sfhomes[order(sfhomes$totvalue, decreasing = FALSE),]
ggplot() +
geom_point(data=sfhomes_low2high, aes(lon,lat, col=totvalue)) +
coord_map()
Try it - Does the output map look different from previous one?
The order of the data in the data frame changes the map display!
Map the sfhomes data in decreasing order by totvalue.
totvaluesfhomes_high2low <- sfhomes[order(sfhomes$totvalue, decreasing = T),]
ggplot() + geom_point(data=sfhomes_high2low, aes(lon,lat, col=totvalue)) +
coord_map()
ggplot GoodnessWhat does this code do?
sfhomes2010_15 <- subset(sfhomes_low2high, as.numeric(SalesYear) > 2009)
ggplot() +
geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 ) +
facet_wrap(~ SalesYear)
ggplot GoodnessVisual spatial analysis!
ggmap extends ggplotCreate basemaps on which you can display your data.
Geocode place names and addresses to get point coordinates.
and more…
Load the libary
library(ggmap)
## Google Maps API Terms of Service: http://developers.google.com/maps/terms.
## Please cite ggmap if you use it: see citation("ggmap") for details.
Some ggmap functionality may require you to register a Google API key
with ggmap
#register_google(key="AIzXXXXXXXXXXXXXXXXXxPSE") # your key here
register_google(key="AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo")
This key will be disabled after the workshop!
ggmapIf you have problems with ggmap you may need to reinstall the library from github.
See this StackOverflow discussion
#devtools::install_github("dkahle/ggmap")
#library(ggmap)
Fetch map data to plot using get_map
Use ?get_map for details
sf_map <- get_map("San Francisco, CA")
This use of the command leverages the Google Geocoding API to determine the coordinates of the named location.
Fetch map data to plot using get_map
Use ?get_map for details
sf_map <- get_map("San Francisco, CA")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=San+Francisco,+CA&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=San%20Francisco%2C%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
ggmap(sf_map)
ggmap(sf_map)
ggmap with point overlaySyntax similar to ggplot
# ggplot() +
ggmap(sf_map) +
geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))
ggmap with point overlayggmap(sf_map) +
geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))
get_mapCreate a basemap zoomed to the extent of our data
# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)
# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes$lon), lat = mean(sfhomes$lat))
sf_ctr # take a look
# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)
get_map ExtentGet a basemap zoomed to the extent of the 2015 data
# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)
# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes15$lon), lat = mean(sfhomes15$lat))
sf_ctr # take a look
## lon lat
## -122.43289 37.75964
# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.759644,-122.432893&zoom=12&size=640x640&scale=1&maptype=terrain&language=en-EN&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
Create a ggmap that uses the customized sf_basemap
sf_basemapggmap(sf_basemap) +
geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))
Finally, let’s add another geospatial data layer to our ggmap.
You can use this method to add as many layers as you want to a ggmap or ggplot.
Use the read.csv function to read in a file of Bart Station locations.
bart <- read.csv("./data/bart.csv")
# take a look
head (bart)
## X Y STATION OPERATOR DIST CO
## 1 -122.2833 37.87406 NORTH BERKELEY BART 4 ALA
## 2 -122.2682 37.86969 DOWNTOWN BERKELEY BART 4 ALA
## 3 -122.2701 37.85321 ASHBY BART 4 ALA
## 4 -122.2518 37.84451 ROCKRIDGE BART 4 ALA
## 5 -122.2671 37.82871 MACARTHUR BART 4 ALA
## 6 -122.2684 37.80808 19TH STREET/OAKLAND BART 4 ALA
ggmap(sf_basemap) +
geom_point(data=sfhomes15, aes(x=lon, y=lat)) +
geom_point(data=bart, aes(x=X,y=Y), col="red")
## Warning: Removed 34 rows containing missing values (geom_point).
Questons?
Determining the geographic coordinates for one or more addresses, zip codes, or named places.
police_stations <- read.csv("data/sf_police_addresses.csv",
stringsAsFactors = F)
head(police_stations)
## PoliceDistrict Phone Address City State
## 1 Central 415-315-2400 766 Vallejo Street San Francisco CA
## 2 Southern 415-575-6000 1251 3rd Street San Francisco CA
## 3 Bayview 415-671-2300 201 Williams Avenue San Francisco CA
## 4 Mission 415-558-5400 630 Valencia Street San Francisco CA
## 5 Northern 415-614-3400 1125 Fillmore Street San Francisco CA
## 6 Park 415-242-3000 1899 Waller Street San Francisco CA
#?geocode
geocode(police_stations)
police_stations$full_address <- paste(police_stations$Address,
police_stations$City, police_stations$State)
head(police_stations)
## PoliceDistrict Phone Address City State
## 1 Central 415-315-2400 766 Vallejo Street San Francisco CA
## 2 Southern 415-575-6000 1251 3rd Street San Francisco CA
## 3 Bayview 415-671-2300 201 Williams Avenue San Francisco CA
## 4 Mission 415-558-5400 630 Valencia Street San Francisco CA
## 5 Northern 415-614-3400 1125 Fillmore Street San Francisco CA
## 6 Park 415-242-3000 1899 Waller Street San Francisco CA
## full_address
## 1 766 Vallejo Street San Francisco CA
## 2 1251 3rd Street San Francisco CA
## 3 201 Williams Avenue San Francisco CA
## 4 630 Valencia Street San Francisco CA
## 5 1125 Fillmore Street San Francisco CA
## 6 1899 Waller Street San Francisco CA
station_coords <- geocode(police_stations$full_address)
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=766%20Vallejo%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1251%203rd%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=201%20Williams%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=630%20Valencia%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1125%20Fillmore%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1899%20Waller%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=461%206th%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1%20Sgt.%20John%20V.%20Young%20Lane%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2345%2024th%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=301%20Eddy%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo
police_stations <- cbind(police_stations, station_coords)
police_stations[1:5,c("PoliceDistrict","full_address","lon","lat")]
## PoliceDistrict full_address lon lat
## 1 Central 766 Vallejo Street San Francisco CA -122.4100 37.79872
## 2 Southern 1251 3rd Street San Francisco CA -122.3894 37.77238
## 3 Bayview 201 Williams Avenue San Francisco CA -122.3980 37.72976
## 4 Mission 630 Valencia Street San Francisco CA -122.4220 37.76285
## 5 Northern 1125 Fillmore Street San Francisco CA -122.4325 37.78022
to your last map
Use a different color
ggmap(sf_basemap) +
geom_point(data=sfhomes15, aes(x=lon, y=lat)) +
geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
geom_point(data=police_stations, aes(x=lon,y=lat), shape=22,
col="black", fill="grey", size=4)
## Warning: Removed 34 rows containing missing values (geom_point).
## About Google geocoding
You can geocode addresses or place names
Super accurate
You may need an API key
Limited to 2500 per day
Use ESRI Geocoding tools if you have more addresses
SF Landmarks
data/landmarks.csv
landmarks <- read.csv("./data/landmarks.csv")
head(landmarks)
## X Y id name
## 1 -13626151 4551548 1 Coit Tower
## 2 -13630804 4544523 2 Twin Peaks
## 3 -13627654 4548289 3 City Hall
## 4 -13635101 4546904 4 Golden Gate Park
ggmap(sf_basemap) +
geom_point(data=sfhomes15, aes(x=lon, y=lat)) +
geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
geom_point(data=landmarks, aes(x=X,y=Y), shape=22,
col="black", fill="grey", size=4)
All good? If not, why not?
## Warning: Removed 34 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
There are limits to what you can do with geospatial data stored in a data frame.
These data are not in geographic coordinates - longitude and latitude.
You can’t map these with ggmap.
## X Y id name
## 1 -13626151 4551548 1 Coit Tower
## 2 -13630804 4544523 2 Twin Peaks
## 3 -13627654 4548289 3 City Hall
## 4 -13635101 4546904 4 Golden Gate Park
The ESRI Shapefile is the most common file format for geospatial data.
Can’t read in or directly plot data in a shapefile with ggplot or ggmap.
Spatial data objects and methods are needed to answer questions like
What properties are in the Noe Valley neighborhood?
What is the average property value in each SF neighborhood?
What is the area of each SF Neighborhood and the property density?
What properties are within walking distance (.25 miles) of the Mission neighborhood?
sp PackageClasses and Methods for Spatial Data
The SP package is most commonly used to construct and manipulate spatial data objects in R.
Hundreds of other R packages that do things with spatial data typically build on SP objects.
sp packageLoad the library sp
Take a look at the different types of spatial object classes supported by sp
library(sp)
getClass("Spatial")
## Warning: package 'sp' was built under R version 3.4.3
## Class "Spatial" [package "sp"]
##
## Slots:
##
## Name: bbox proj4string
## Class: matrix CRS
##
## Known Subclasses:
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2
| Geometry | Spatial Object | Spatial Object with Attributes |
|---|---|---|
| Points | SpatialPoints | SpatialPointsDataFrame |
| Lines | SpatialLines | SpatialLinesDataFrame |
| Polygons | SpatialPolygons | SpatialPolygonsDataFrame |
We use the S*DF objects most frequently!
Let’s transform the sfhomes15 data frame to an sp object of type SpatialPointsDataFrame
sp::coordinates()Use the sp::coordinates() method
When transforming a DF to SPDF you are setting the coordinates.
Note - the transformation happens in place!
sfhomes15_sp <- sfhomes15 # Make a copy - why?
class(sfhomes15_sp) # check the class of the object
coordinates(sfhomes15_sp) <- c('lon','lat') # Make it spatial - ORDER MATTERS!!
class(sfhomes15_sp) # check the class of the object
sfhomes15_sp <- sfhomes15 # Make a copy - why?
class(sfhomes15_sp) # check the class of the object
## [1] "data.frame"
coordinates(sfhomes15_sp) <- c('lon','lat') # Make it spatial - ORDER MATTERS!!
class(sfhomes15_sp) # check the class of the object
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
You transformed a data frame to an SPDF using:
coordinates(sfhomes15_sp) <- c('lon','lat')
Now try this:
coordinates(sfhomes15_sp)
coordinates(sfhomes15_sp)
## lon lat
## 24 -122.4741 37.73601
## 35 -122.3965 37.75920
## 76 -122.3942 37.77181
## 79 -122.4754 37.71518
## 85 -122.4210 37.75876
## 90 -122.4149 37.79362
## 92 -122.4370 37.74020
## 98 -122.4437 37.78881
## 223 -122.4754 37.71492
## 229 -122.3919 37.78868
## 273 -122.3968 37.77945
## 341 -122.5040 37.75515
## 363 -122.4041 37.76422
## 366 -122.3934 37.77764
## 388 -122.4166 37.74931
## 451 -122.4112 37.73683
## 467 -122.4507 37.76448
## 473 -122.3990 37.75450
## 487 -122.4195 37.76649
## 534 -122.4225 37.78676
## 548 -122.4410 37.70854
## 549 -122.5009 37.74485
## 550 -122.4220 37.73048
## 559 -122.4947 37.75832
## 567 -122.3942 37.77181
## 595 -122.4187 37.80469
## 624 -122.4019 37.77323
## 631 -122.4229 37.80586
## 646 -122.4225 37.78676
## 658 -122.3919 37.78868
## 665 -122.4261 37.74349
## 668 -122.3904 37.78139
## 683 -122.4458 37.74398
## 774 -122.3883 37.73754
## 776 -122.4767 37.75347
## 864 -122.4326 37.74229
## 974 -122.4380 37.72624
## 1018 -122.4452 37.76906
## 1097 -122.4937 37.77999
## 1155 -122.5047 37.75390
## 1164 -122.4498 37.75451
## 1167 -122.4055 37.80424
## 1188 -122.4194 37.75596
## 1242 -122.4088 37.77953
## 1273 -122.4331 37.76019
## 1374 -122.3918 37.77899
## 1379 -122.4070 37.75674
## 1487 -122.4236 37.77800
## 1598 -122.4758 37.75497
## 1600 -122.4419 37.74582
## 1618 -122.5030 37.74587
## 1668 -122.3989 37.78990
## 1698 -122.4236 37.77800
## 1747 -122.3972 37.76314
## 1755 -122.4434 37.75159
## 1768 -122.4242 37.80551
## 1804 -122.4549 37.76402
## 1828 -122.4320 37.78195
## 1833 -122.4405 37.71825
## 1913 -122.4146 37.74714
## 1999 -122.3942 37.77181
## 2024 -122.4676 37.71048
## 2033 -122.4249 37.76548
## 2066 -122.4243 37.75206
## 2073 -122.4097 37.77613
## 2096 -122.3996 37.73859
## 2248 -122.4045 37.79929
## 2258 -122.4442 37.78322
## 2285 -122.3919 37.78868
## 2314 -122.4026 37.77879
## 2329 -122.4368 37.78228
## 2339 -122.4240 37.74366
## 2355 -122.4126 37.79035
## 2416 -122.4579 37.72632
## 2428 -122.4089 37.79286
## 2436 -122.4409 37.80451
## 2443 -122.4127 37.78995
## 2500 -122.4357 37.76468
## 2508 -122.4666 37.76347
## 2537 -122.4870 37.74605
## 2538 -122.3942 37.77181
## 2627 -122.4219 37.77389
## 2692 -122.3958 37.77705
## 2717 -122.4116 37.77806
## 2734 -122.4290 37.72977
## 2793 -122.4764 37.73390
## 2914 -122.4183 37.76712
## 2969 -122.4711 37.72241
## 2975 -122.4393 37.78018
## 3127 -122.4393 37.78951
## 3157 -122.3702 37.72838
## 3205 -122.4398 37.78746
## 3262 -122.4540 37.78251
## 3288 -122.4698 37.75298
## 3296 -122.4607 37.73104
## 3306 -122.4081 37.77962
## 3354 -122.4849 37.74663
## 3403 -122.4209 37.75792
## 3436 -122.4854 37.74788
## 3471 -122.4856 37.78116
## 3485 -122.4177 37.80345
## 3516 -122.4354 37.72201
## 3578 -122.4199 37.77606
## 3611 -122.4269 37.80229
## 3658 -122.4108 37.76534
## 3754 -122.3888 37.72872
## 3772 -122.3891 37.73264
## 4023 -122.4827 37.77762
## 4025 -122.4450 37.76430
## 4058 -122.4328 37.71546
## 4076 -122.4660 37.74828
## 4111 -122.5023 37.74333
## 4114 -122.4970 37.74468
## 4140 -122.4741 37.73601
## 4150 -122.4060 37.79853
## 4249 -122.4302 37.74726
## 4270 -122.4194 37.76710
## 4287 -122.4676 37.71048
## 4369 -122.4950 37.73327
## 4376 -122.3919 37.78868
## 4393 -122.3972 37.77570
## 4427 -122.4332 37.77869
## 4466 -122.4786 37.77394
## 4469 -122.4210 37.76470
## 4483 -122.4318 37.71761
## 4514 -122.4214 37.78139
## 4526 -122.4153 37.74047
## 4527 -122.4209 37.79368
## 4555 -122.4145 37.79167
## 4561 -122.3910 37.78938
## 4566 -122.4031 37.79058
## 4578 -122.3906 37.75501
## 4579 -122.4154 37.74537
## 4592 -122.3922 37.78607
## 4599 -122.4385 37.75980
## 4602 -122.4785 37.76295
## 4618 -122.4137 37.79099
## 4648 -122.4635 37.76461
## 4663 -122.4160 37.74563
## 4669 -122.3918 37.77899
## 4759 -122.4911 37.75318
## 4800 -122.4137 37.79099
## 4825 -122.4147 37.72528
## 4853 -122.4388 37.75978
## 4864 -122.3696 37.72900
## 4896 -122.4245 37.79493
## 4899 -122.3902 37.75236
## 4921 -122.5103 37.77369
## 4939 -122.3930 37.73299
## 4944 -122.3919 37.78868
## 4970 -122.4646 37.77383
## 4986 -122.3904 37.78139
## 5053 -122.4287 37.77484
## 5085 -122.4244 37.78587
## 5104 -122.4255 37.72744
## 5108 -122.4364 37.77781
## 5109 -122.4394 37.80326
## 5114 -122.4853 37.75287
## 5188 -122.4302 37.74806
## 5228 -122.3881 37.73174
## 5247 -122.4242 37.80551
## 5252 -122.4296 37.79988
## 5261 -122.4048 37.79879
## 5287 -122.4335 37.79946
## 5288 -122.4611 37.73481
## 5316 -122.3918 37.77899
## 5320 -122.4481 37.77945
## 5321 -122.4477 37.71766
## 5334 -122.3919 37.78868
## 5350 -122.4336 37.72535
## 5351 -122.4540 37.78251
## 5360 -122.4264 37.72624
## 5366 -122.3922 37.78607
## 5398 -122.3905 37.78759
## 5425 -122.3919 37.78868
## 5430 -122.4263 37.72116
## 5450 -122.4175 37.79865
## 5461 -122.3989 37.78577
## 5467 -122.4175 37.77341
## 5511 -122.3919 37.78868
## 5519 -122.4093 37.74024
## 5521 -122.4816 37.77294
## 5522 -122.4194 37.75596
## 5547 -122.3959 37.79054
## 5561 -122.4213 37.78884
## 5571 -122.4194 37.76710
## 5609 -122.4245 37.79920
## 5647 -122.4085 37.75911
## 5654 -122.4673 37.75439
## 5671 -122.4459 37.75811
## 5691 -122.4806 37.77724
## 5705 -122.4167 37.76091
## 5712 -122.3702 37.72838
## 5732 -122.4207 37.73596
## 5740 -122.4101 37.77808
## 5757 -122.4752 37.71513
## 5765 -122.3922 37.73203
## 5827 -122.4036 37.71724
## 5978 -122.4264 37.73053
## 5997 -122.3886 37.77159
## 6043 -122.5068 37.77785
## 6056 -122.4089 37.79286
## 6134 -122.4191 37.70898
## 6140 -122.4107 37.76483
## 6181 -122.3974 37.78575
## 6191 -122.4321 37.73585
## 6208 -122.4774 37.75583
## 6252 -122.4355 37.74743
## 6258 -122.3904 37.78271
## 6342 -122.4846 37.75914
## 6365 -122.4356 37.75016
## 6410 -122.3824 37.72887
## 6430 -122.4699 37.75813
## 6477 -122.4656 37.71671
## 6581 -122.4792 37.77343
## 6687 -122.4243 37.79352
## 6691 -122.4371 37.71173
## 6694 -122.4399 37.71837
## 6721 -122.4308 37.73647
## 6723 -122.4566 37.75242
## 6782 -122.4298 37.71104
## 6787 -122.4405 37.80504
## 6814 -122.4848 37.73940
## 6863 -122.4122 37.79293
## 6892 -122.4419 37.74582
## 6903 -122.4084 37.79154
## 6915 -122.4144 37.77516
## 6918 -122.3958 37.77816
## 6932 -122.4244 37.74805
## 6933 -122.4368 37.78249
## 6988 -122.5064 37.77555
## 7006 -122.4528 37.76708
## 7009 -122.4947 37.76067
## 7032 -122.4050 37.77762
## 7073 -122.4194 37.75596
## 7142 -122.4953 37.73096
## 7149 -122.4104 37.75357
## 7175 -122.4370 37.76934
## 7206 -122.4936 37.77938
## 7303 -122.4327 37.77299
## 7320 -122.4194 37.76710
## 7337 -122.4752 37.71551
## 7395 -122.4263 37.76856
## 7404 -122.3939 37.75915
## 7447 -122.4113 37.72850
## 7462 -122.4127 37.76517
## 7493 -122.4194 37.75596
## 7533 -122.4515 37.75418
## 7555 -122.4071 37.77553
## 7564 -122.4389 37.73129
## 7600 -122.4097 37.77613
## 7649 -122.4290 37.78872
## 7651 -122.3958 37.77816
## 7704 -122.4337 37.76214
## 7735 -122.3904 37.78271
## 7832 -122.4853 37.75355
## 7833 -122.3983 37.75753
## 7848 -122.4490 37.76795
## 7936 -122.4087 37.77310
## 7981 -122.4489 37.71967
## 7991 -122.4325 37.77805
## 8047 -122.4019 37.78005
## 8089 -122.4194 37.75596
## 8127 -122.4685 37.76188
## 8161 -122.4307 37.74163
## 8275 -122.4126 37.77870
## 8322 -122.3924 37.73700
## 8344 -122.4126 37.77870
## 8378 -122.4900 37.78324
## 8401 -122.4642 37.71638
## 8410 -122.4279 37.77408
## 8419 -122.4016 37.75145
## 8422 -122.4836 37.71006
## 8438 -122.4668 37.73426
## 8517 -122.3919 37.78868
## 8537 -122.3980 37.71952
## 8590 -122.4291 37.73289
## 8612 -122.4194 37.75596
## 8652 -122.4222 37.72570
## 8694 -122.3919 37.78868
## 8703 -122.3958 37.77339
## 8713 -122.4301 37.77286
## 8730 -122.4126 37.77870
## 8818 -122.3897 37.78278
## 8823 -122.4361 37.78937
## 8842 -122.4257 37.79707
## 8844 -122.4259 37.78673
## 8888 -122.4318 37.79685
## 8891 -122.4277 37.77402
## 8899 -122.3896 37.78063
## 8907 -122.4728 37.75198
## 8932 -122.4102 37.74514
## 8962 -122.3942 37.77181
## 8983 -122.4206 37.80236
## 9049 -122.4651 37.75119
## 9054 -122.4371 37.73361
## 9179 -122.4262 37.72086
## 9216 -122.4160 37.73419
## 9219 -122.4550 37.74285
## 9356 -122.4522 37.73356
## 9376 -122.3982 37.75309
## 9384 -122.4152 37.73405
## 9469 -122.4174 37.80442
## 9479 -122.4251 37.72966
## 9490 -122.4835 37.74208
## 9569 -122.4723 37.73269
## 9639 -122.4883 37.77411
## 9640 -122.4509 37.73126
## 9647 -122.3909 37.75967
## 9675 -122.3891 37.76366
## 9691 -122.4381 37.78210
## 9764 -122.4875 37.76255
## 9776 -122.4422 37.72376
## 9806 -122.3962 37.77393
## 9830 -122.4642 37.71673
## 9906 -122.4433 37.76135
## 9959 -122.3934 37.77764
## 9966 -122.4639 37.72971
## 9985 -122.4411 37.75920
## 10010 -122.4279 37.76639
## 10046 -122.4755 37.71495
## 10063 -122.3879 37.73919
## 10066 -122.4429 37.78894
## 10100 -122.3890 37.73260
## 10157 -122.4377 37.79704
## 10164 -122.4180 37.74667
## 10170 -122.4453 37.79137
## 10229 -122.4137 37.79099
## 10239 -122.4451 37.78484
## 10257 -122.4102 37.74590
## 10263 -122.4353 37.75147
## 10272 -122.4974 37.78118
## 10273 -122.4401 37.75541
## 10307 -122.4167 37.79138
## 10338 -122.4343 37.73743
## 10343 -122.4108 37.77597
## 10345 -122.5037 37.74846
## 10351 -122.3910 37.78938
## 10353 -122.4865 37.76343
## 10410 -122.3940 37.75300
## 10415 -122.4287 37.73629
## 10419 -122.4344 37.76357
## 10499 -122.4531 37.73537
## 10530 -122.4150 37.73959
## 10578 -122.4594 37.73954
## 10583 -122.4760 37.71558
## 10600 -122.4832 37.78375
## 10607 -122.4911 37.73160
## 10632 -122.4347 37.78713
## 10646 -122.4518 37.76347
## 10661 -122.4066 37.73036
## 10781 -122.4124 37.80434
## 10782 -122.4554 37.71773
## 10800 -122.4152 37.73705
## 10857 -122.4178 37.76399
## 10887 -122.4260 37.74559
## 10914 -122.4349 37.74579
## 10957 -122.4173 37.74421
## 10967 -122.4587 37.73858
## 10985 -122.4325 37.77540
## 11005 -122.4558 37.73265
## 11023 -122.4173 37.73753
## 11035 -122.4386 37.76810
## 11094 -122.4275 37.72671
## 11158 -122.3937 37.78702
## 11160 -122.4236 37.77800
## 11170 -122.5054 37.76363
## 11171 -122.4001 37.79803
## 11180 -122.4757 37.71556
## 11181 -122.4359 37.75531
## 11202 -122.4133 37.78070
## 11258 -122.4437 37.80053
## 11265 -122.4865 37.76465
## 11284 -122.4117 37.75704
## 11286 -122.4194 37.75596
## 11305 -122.4136 37.79298
## 11385 -122.4424 37.74562
## 11451 -122.4370 37.73761
## 11498 -122.4438 37.71852
## 11554 -122.4990 37.77913
## 11564 -122.4194 37.75596
## 11576 -122.3910 37.78938
## 11598 -122.3824 37.72852
## 11603 -122.4748 37.71490
## 11622 -122.3988 37.71372
## 11636 -122.4678 37.72225
## 11690 -122.4415 37.78044
## 11728 -122.4793 37.75869
## 11729 -122.4165 37.73373
## 11808 -122.4457 37.75769
## 11830 -122.4965 37.77600
## 11844 -122.4468 37.73736
## 11873 -122.4474 37.73477
## 11876 -122.4126 37.79035
## 11882 -122.3970 37.75975
## 11918 -122.4438 37.76155
## 11927 -122.4333 37.73432
## 11968 -122.3936 37.78370
## 11969 -122.3907 37.78117
## 11971 -122.4742 37.75636
## 11972 -122.4542 37.77443
## 12052 -122.4385 37.75442
## 12232 -122.5060 37.73665
## 12285 -122.4786 37.77955
## 12350 -122.4672 37.73270
## 12383 -122.4093 37.77680
## 12428 -122.4274 37.77132
## 12467 -122.3971 37.78091
## 12533 -122.4024 37.77546
## 12559 -122.4130 37.80397
## 12569 -122.4114 37.72936
## 12622 -122.5059 37.73682
## 12668 -122.4131 37.79777
## 12672 -122.4154 37.73877
## 12675 -122.4357 37.73740
## 12730 -122.4011 37.73263
## 12770 -122.4221 37.72534
## 12783 -122.3910 37.78938
## 12865 -122.3919 37.78868
## 12866 -122.4891 37.75507
## 13000 -122.3968 37.77945
## 13015 -122.4341 37.76397
## 13054 -122.4165 37.77528
## 13177 -122.4043 37.78622
## 13251 -122.4294 37.77067
## 13252 -122.3933 37.77891
## 13268 -122.4316 37.76576
## 13293 -122.4936 37.77938
## 13324 -122.4116 37.75238
## 13382 -122.4513 37.76809
## 13391 -122.4248 37.80057
## 13414 -122.4433 37.78314
## 13422 -122.4141 37.80393
## 13470 -122.4439 37.71527
## 13472 -122.4815 37.74502
## 13478 -122.4214 37.78139
## 13507 -122.4298 37.73805
## 13529 -122.4242 37.80551
## 13550 -122.4439 37.80384
## 13581 -122.3956 37.73468
## 13593 -122.4286 37.74623
## 13595 -122.4617 37.78577
## 13662 -122.3918 37.77899
## 13675 -122.3696 37.72900
## 13701 -122.3962 37.77393
## 13722 -122.4025 37.73527
## 13777 -122.4750 37.71499
## 13821 -122.4032 37.78824
## 13895 -122.4485 37.73998
## 13903 -122.3919 37.78868
## 13926 -122.4141 37.75915
## 13936 -122.4296 37.73278
## 14044 -122.4484 37.76533
## 14072 -122.4546 37.74500
## 14087 -122.4642 37.75929
## 14239 -122.4665 37.71392
## 14241 -122.4456 37.80221
## 14256 -122.4990 37.77356
## 14287 -122.4234 37.77253
## 14338 -122.4475 37.77957
## 14345 -122.4192 37.75400
## 14353 -122.3696 37.72900
## 14367 -122.4253 37.79247
## 14413 -122.4397 37.76158
## 14496 -122.4276 37.79552
## 14517 -122.4617 37.75037
## 14679 -122.4565 37.72132
## 14690 -122.4365 37.71167
## 14699 -122.4375 37.76305
## 14721 -122.3919 37.78868
## 14728 -122.3934 37.77764
## 14769 -122.4258 37.78636
## 14773 -122.4326 37.75640
## 14777 -122.3936 37.78370
## 14803 -122.4121 37.77255
## 14870 -122.4481 37.73958
## 14965 -122.4093 37.74690
## 14974 -122.5070 37.78027
## 15069 -122.4340 37.74511
## 15076 -122.4733 37.77618
## 15143 -122.4234 37.77253
## 15159 -122.4937 37.78119
## 15167 -122.4129 37.78791
## 15179 -122.4668 37.74374
## 15186 -122.4476 37.77794
## 15187 -122.3991 37.78318
## 15197 -122.3943 37.77898
## 15208 -122.4557 37.71642
## 15217 -122.4678 37.72282
## 15260 -122.4416 37.71429
## 15267 -122.4507 37.74389
## 15319 -122.4345 37.78565
## 15327 -122.4957 37.77863
## 15328 -122.4564 37.70934
## 15349 -122.4151 37.73756
## 15357 -122.3761 37.73193
## 15426 -122.4009 37.76461
## 15474 -122.4009 37.76461
## 15485 -122.5046 37.75155
## 15486 -122.4718 37.76104
## 15498 -122.4174 37.74938
## 15525 -122.4124 37.77093
## 15551 -122.3943 37.78613
## 15567 -122.4055 37.80424
## 15580 -122.3914 37.78451
## 15587 -122.5046 37.74607
## 15589 -122.4339 37.72119
## 15604 -122.4248 37.75359
## 15615 -122.4098 37.71496
## 15621 -122.4667 37.73841
## 15666 -122.4379 37.75930
## 15669 -122.4194 37.75596
## 15713 -122.3702 37.72838
## 15717 -122.4382 37.75646
## 15790 -122.4021 37.76545
## 15887 -122.4311 37.76352
## 15918 -122.3942 37.77181
## 15933 -122.4178 37.76399
## 15945 -122.4390 37.79617
## 16033 -122.4014 37.78629
## 16076 -122.4149 37.79362
## 16105 -122.4489 37.76179
## 16173 -122.4135 37.72849
## 16199 -122.4014 37.77974
## 16204 -122.4327 37.76614
## 16210 -122.4210 37.75876
## 16212 -122.4126 37.77870
## 16215 -122.3958 37.77705
## 16229 -122.4400 37.71831
## 16310 -122.3988 37.72731
## 16315 -122.4254 37.80339
## 16389 -122.4800 37.75339
## 16440 -122.4970 37.74495
## 16450 -122.4194 37.75596
## 16458 -122.4624 37.77409
## 16474 -122.4542 37.77443
## 16480 -122.4014 37.77974
## 16500 -122.4126 37.77870
## 16507 -122.4300 37.78804
## 16517 -122.4240 37.76628
## 16531 -122.4266 37.73311
## 16549 -122.5038 37.74097
## 16640 -122.4384 37.77517
## 16643 -122.4517 37.74504
## 16712 -122.4194 37.75596
## 16755 -122.3895 37.73005
## 16778 -122.4708 37.75810
## 16796 -122.4035 37.72002
## 16800 -122.3696 37.72900
## 16802 -122.4042 37.71115
## 16921 -122.3696 37.72900
## 16934 -122.4054 37.77794
## 17006 -122.4459 37.77465
## 17051 -122.4560 37.75650
## 17073 -122.4386 37.73261
## 17084 -122.3942 37.77181
## 17091 -122.4750 37.71499
## 17156 -122.4286 37.71810
## 17168 -122.4766 37.76075
## 17199 -122.3910 37.78938
## 17200 -122.4151 37.78524
## 17223 -122.4011 37.77047
## 17249 -122.4159 37.74769
## 17252 -122.4908 37.73143
## 17265 -122.4043 37.80443
## 17304 -122.4014 37.77974
## 17345 -122.4245 37.79493
## 17358 -122.4704 37.72549
## 17384 -122.3904 37.78271
## 17404 -122.4080 37.77576
## 17410 -122.3904 37.75368
## 17465 -122.5033 37.74841
## 17491 -122.4015 37.73255
## 17510 -122.4359 37.75497
## 17517 -122.4519 37.76025
## 17542 -122.4124 37.77093
## 17562 -122.4514 37.70881
## 17613 -122.3893 37.73988
## 17618 -122.4249 37.75375
## 17698 -122.4262 37.79417
## 17715 -122.4140 37.75590
## 17754 -122.4360 37.80400
## 17767 -122.4916 37.78122
## 17857 -122.4065 37.77812
## 17867 -122.4777 37.76001
## 17876 -122.5107 37.77399
## 17881 -122.4226 37.73575
## 17887 -122.4188 37.77159
## 17912 -122.3886 37.77159
## 17929 -122.4289 37.72161
## 17933 -122.4452 37.71900
## 17936 -122.4458 37.72772
## 17986 -122.4429 37.73288
## 17988 -122.4473 37.73331
## 18002 -122.4422 37.77803
## 18018 -122.4772 37.74409
## 18050 -122.4198 37.76161
## 18053 -122.4038 37.76405
## 18059 -122.4221 37.74746
## 18071 -122.4271 37.73809
## 18107 -122.4388 37.73398
## 18164 -122.4160 37.73504
## 18227 -122.4805 37.75531
## 18231 -122.4671 37.76370
## 18323 -122.4229 37.80586
## 18363 -122.3727 37.73057
## 18386 -122.4767 37.77714
## 18433 -122.4408 37.78769
## 18439 -122.4345 37.73587
## 18444 -122.4716 37.78657
## 18453 -122.3928 37.78205
## 18470 -122.4682 37.75544
## 18504 -122.5075 37.74854
## 18511 -122.3959 37.79054
## 18553 -122.3942 37.77181
## 18584 -122.4140 37.75590
## 18619 -122.4315 37.73498
## 18627 -122.4716 37.75516
## 18642 -122.4780 37.78258
## 18653 -122.4258 37.72884
## 18718 -122.4391 37.75547
## 18733 -122.4294 37.72902
## 18742 -122.4194 37.75596
## 18792 -122.4181 37.78824
## 18795 -122.3938 37.78638
## 18801 -122.4288 37.72584
## 18871 -122.4234 37.77253
## 18886 -122.4305 37.74402
## 18913 -122.4280 37.74869
## 18960 -122.4387 37.71054
## 18968 -122.4041 37.80347
## 18971 -122.4058 37.80159
## 18990 -122.4398 37.75879
## 19015 -122.4613 37.76457
## 19065 -122.4121 37.77187
## 19082 -122.4684 37.76042
## 19105 -122.4836 37.71006
## 19143 -122.4642 37.78458
## 19153 -122.4520 37.73468
## 19197 -122.4102 37.75679
## 19202 -122.4412 37.71904
## 19227 -122.3896 37.78063
## 19281 -122.4309 37.72308
## 19325 -122.4287 37.77642
## 19375 -122.4152 37.73366
## 19532 -122.4596 37.78655
## 19535 -122.5019 37.77481
## 19570 -122.4623 37.73544
## 19626 -122.4348 37.74426
## 19641 -122.3919 37.78868
## 19718 -122.4194 37.75596
## 19765 -122.4194 37.75596
## 19782 -122.3911 37.73912
## 19802 -122.4152 37.79845
## 19811 -122.3951 37.72508
## 19832 -122.4798 37.78437
## 19838 -122.4301 37.74120
## 19845 -122.4435 37.75964
## 19880 -122.4243 37.79352
## 19897 -122.3934 37.77764
## 19932 -122.4590 37.71972
## 19933 -122.4261 37.73585
## 19954 -122.4934 37.76073
## 19991 -122.4242 37.79788
## 20009 -122.4528 37.71555
## 20026 -122.3909 37.75967
## 20061 -122.4382 37.80320
## 20095 -122.3919 37.78868
## 20102 -122.4936 37.74282
## 20221 -122.4251 37.75560
## 20235 -122.4709 37.75704
## 20293 -122.4937 37.77999
## 20310 -122.4252 37.74835
## 20320 -122.4355 37.71758
## 20321 -122.4739 37.74485
## 20347 -122.5013 37.76014
## 20356 -122.4517 37.73356
## 20426 -122.4014 37.77974
## 20484 -122.4345 37.80458
## 20562 -122.3974 37.78575
## 20605 -122.4729 37.77471
## 20611 -122.3922 37.78607
## 20655 -122.4621 37.77642
## 20687 -122.4548 37.78306
## 20743 -122.4469 37.71901
## 20754 -122.4493 37.77461
## 20768 -122.3962 37.77393
## 20778 -122.4743 37.76473
## 20779 -122.4024 37.76117
## 20804 -122.5022 37.74887
## 20871 -122.3922 37.73197
## 20877 -122.4180 37.76654
## 20928 -122.4731 37.77999
## 20971 -122.4816 37.77795
## 20979 -122.4582 37.73860
## 21007 -122.3696 37.72900
## 21053 -122.4194 37.75596
## 21064 -122.4263 37.76856
## 21079 -122.4425 37.80069
## 21093 -122.4336 37.74361
## 21105 -122.4245 37.79920
## 21124 -122.4399 37.79230
## 21152 -122.4642 37.72043
## 21174 -122.4227 37.79931
## 21201 -122.4340 37.77121
## 21231 -122.4272 37.72286
## 21255 -122.4094 37.75902
## 21374 -122.4208 37.73601
## 21424 -122.3702 37.72838
## 21455 -122.4195 37.80612
## 21456 -122.3932 37.78006
## 21458 -122.4980 37.75912
## 21479 -122.4331 37.79233
## 21486 -122.3922 37.78607
## 21496 -122.4739 37.78414
## 21499 -122.4031 37.79058
## 21508 -122.4606 37.72592
## 21522 -122.4188 37.77159
## 21599 -122.4630 37.74288
## 21626 -122.4271 37.78448
## 21632 -122.4368 37.74525
## 21644 -122.4176 37.76271
## 21652 -122.3901 37.75370
## 21678 -122.3702 37.72838
## 21694 -122.4414 37.79054
## 21711 -122.4126 37.77870
## 21712 -122.3904 37.78139
## 21790 -122.4446 37.79053
## 21791 -122.5024 37.73799
## 21834 -122.4218 37.79478
## 21889 -122.4157 37.79792
## 21932 -122.4520 37.71732
## 21957 -122.4282 37.79869
## 21972 -122.4552 37.74086
## 21985 -122.4549 37.73345
## 22007 -122.4195 37.74526
## 22024 -122.5006 37.75075
## 22064 -122.4569 37.75564
## 22085 -122.4497 37.72884
## 22111 -122.3933 37.77891
## 22167 -122.5041 37.73937
## 22254 -122.3992 37.79825
## 22255 -122.4107 37.77166
## 22273 -122.4448 37.71416
## 22324 -122.3986 37.71395
## 22340 -122.4376 37.75052
## 22499 -122.4438 37.72965
## 22531 -122.4141 37.77608
## 22533 -122.4409 37.74887
## 22549 -122.4431 37.71275
## 22554 -122.4009 37.76461
## 22600 -122.4362 37.74701
## 22607 -122.4753 37.71516
## 22673 -122.4131 37.74787
## 22739 -122.3919 37.78868
## 22801 -122.4709 37.75504
## 22802 -122.4014 37.77974
## 22820 -122.4887 37.77819
## 22870 -122.4367 37.76372
## 22871 -122.4507 37.72158
## 22896 -122.4342 37.74955
## 22924 -122.4460 37.71426
## 22931 -122.4210 37.74497
## 22951 -122.4948 37.77410
## 23000 -122.4254 37.76334
## 23036 -122.4163 37.80401
## 23045 -122.3996 37.77438
## 23056 -122.4753 37.71487
## 23074 -122.3696 37.72900
## 23117 -122.4223 37.79267
## 23135 -122.3942 37.77181
## 23154 -122.4648 37.75510
## 23172 -122.4415 37.78044
## 23173 -122.4835 37.78257
## 23251 -122.4119 37.71432
## 23267 -122.4263 37.71022
## 23308 -122.4769 37.76468
## 23356 -122.4613 37.71719
## 23400 -122.4689 37.77570
## 23433 -122.4753 37.71553
## 23467 -122.3928 37.78205
## 23478 -122.4855 37.75523
## 23486 -122.4276 37.72233
## 23503 -122.4108 37.76534
## 23519 -122.4454 37.76390
## 23535 -122.4639 37.76393
## 23540 -122.4599 37.72230
## 23542 -122.4860 37.78680
## 23598 -122.3910 37.78938
## 23601 -122.3922 37.78607
## 23721 -122.4452 37.71600
## 23782 -122.4012 37.75378
## 23790 -122.5038 37.75594
## 23828 -122.3957 37.73636
## 23881 -122.3906 37.75501
## 23895 -122.4517 37.72883
## 23896 -122.4011 37.77047
## 23913 -122.4734 37.75749
## 23918 -122.4676 37.71048
## 23919 -122.4936 37.77938
## 23991 -122.4418 37.73477
## 24130 -122.4544 37.76267
## 24185 -122.4891 37.77817
## 24223 -122.4194 37.75596
## 24282 -122.4934 37.73923
## 24314 -122.4577 37.74434
## 24326 -122.3880 37.73296
## 24374 -122.3696 37.72900
## 24393 -122.4372 37.75934
## 24401 -122.4517 37.73329
## 24402 -122.4711 37.74663
## 24504 -122.4536 37.70915
## 24546 -122.4544 37.78545
## 24554 -122.4395 37.75941
## 24574 -122.3997 37.75216
## 24607 -122.4836 37.77285
## 24639 -122.4824 37.75738
## 24659 -122.3920 37.72865
## 24674 -122.4760 37.78306
## 24700 -122.4138 37.73519
## 24716 -122.4922 37.77477
## 24724 -122.4374 37.75893
## 24742 -122.4704 37.78430
## 24781 -122.4137 37.79099
## 24794 -122.4016 37.75125
## 24816 -122.4009 37.72799
## 24840 -122.4704 37.74936
## 24872 -122.4236 37.77800
## 24896 -122.4441 37.73437
## 24922 -122.4446 37.77193
## 24940 -122.4192 37.72509
## 24949 -122.4356 37.72596
## 24966 -122.3919 37.78868
## 24983 -122.4458 37.73512
## 24986 -122.4668 37.77558
str(sfhomes15) # the data frame
str(sfhomes15_sp) # the data frame
str(sfhomes15) # the data frame
## 'data.frame': 835 obs. of 19 variables:
## $ FiscalYear : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## $ SalesDate : chr "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
## $ Address : chr "0000 2760 19TH AV0015" "0000 0560AMISSOURI ST0000" "0000 0718 LONG BRIDGE ST1202" "0000 0103 SUMMIT WY0000" ...
## $ YearBuilt : int 1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
## $ NumBedrooms : int 2 2 2 3 3 3 3 2 0 2 ...
## $ NumBathrooms : int 2 2 2 3 3 3 3 2 0 2 ...
## $ NumRooms : int 5 5 5 6 7 7 7 5 0 0 ...
## $ NumStories : int 0 1 1 0 1 14 2 1 0 1 ...
## $ NumUnits : int 1 1 1 1 1 1 1 1 0 1 ...
## $ AreaSquareFeet : int 1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
## $ ImprovementValue: int 432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
## $ LandValue : int 432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
## $ Neighborhood : chr "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
## $ Location : chr "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
## $ SupeDistrict : int 7 10 6 7 9 3 8 2 7 6 ...
## $ totvalue : int 865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
## $ SalesYear : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## $ lat : num 37.7 37.8 37.8 37.7 37.8 ...
## $ lon : num -122 -122 -122 -122 -122 ...
str(sfhomes15_sp) # the SPDF
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 835 obs. of 17 variables:
## .. ..$ FiscalYear : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## .. ..$ SalesDate : chr [1:835] "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
## .. ..$ Address : chr [1:835] "0000 2760 19TH AV0015" "0000 0560AMISSOURI ST0000" "0000 0718 LONG BRIDGE ST1202" "0000 0103 SUMMIT WY0000" ...
## .. ..$ YearBuilt : int [1:835] 1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
## .. ..$ NumBedrooms : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
## .. ..$ NumBathrooms : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
## .. ..$ NumRooms : int [1:835] 5 5 5 6 7 7 7 5 0 0 ...
## .. ..$ NumStories : int [1:835] 0 1 1 0 1 14 2 1 0 1 ...
## .. ..$ NumUnits : int [1:835] 1 1 1 1 1 1 1 1 0 1 ...
## .. ..$ AreaSquareFeet : int [1:835] 1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
## .. ..$ ImprovementValue: int [1:835] 432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
## .. ..$ LandValue : int [1:835] 432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
## .. ..$ Neighborhood : chr [1:835] "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
## .. ..$ Location : chr [1:835] "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
## .. ..$ SupeDistrict : int [1:835] 7 10 6 7 9 3 8 2 7 6 ...
## .. ..$ totvalue : int [1:835] 865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
## .. ..$ SalesYear : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## ..@ coords.nrs : int [1:2] 19 18
## ..@ coords : num [1:835, 1:2] -122 -122 -122 -122 -122 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:835] "24" "35" "76" "79" ...
## .. .. ..$ : chr [1:2] "lon" "lat"
## ..@ bbox : num [1:2, 1:2] -122.5 37.7 -122.4 37.8
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "lon" "lat"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NA
You can see from str(sfhomes) that a SPDF object is a collection of slots or components. The key ones are:
@data data frame of attributes that describe each location@coords the coordinates for each geometric object - here points@bbox the min and max bounding coordinates@proj4string the coordinate reference system defintion as a stringExplore the object in the Environment window
Review the output of each of these:
summary(sfhomes15_sp)
head(sfhomes15_sp@data)
class(sfhomes15_sp@data)
sfhomes15_sp@bbox
bbox(sfhomes15_sp)
head(sfhomes15_sp@coords)
head(sfhomes15_sp$lat)
head(sfhomes15_sp$lon)
sfhomes15_sp@proj4string
proj4string(sfhomes15_sp)
Look at sfhomes15_sp in the environment window
Are all the columns that were present in sfhomes15 also in sfhomes15_sp?
Is there a slot in sfhomes15_sp without data?
proj4string(sfhomes15_sp) # get a CRS object
## [1] NA
SpatialPointsDataFrameplot(sfhomes15_sp) # using sp::plot
We created great maps of sfhomes point data with ggplot and ggmap.
Then we created a simple map of the SPDF with plot.
We aren’t seeing the value of sp objects just yet.
rgdalrgdal is an R port of the powerful and widely used GDAL library.
It is the most commonly used R library for importing and exporting spatial data.
OGR: for vector data: readOGR() and writeOGR()
GDAL for raster data: readGDAL() and writeGDAL()
rgdallibrary(rgdal)
## Warning: package 'rgdal' was built under R version 3.4.3
## rgdal: version: 1.2-16, (SVN revision 701)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
## Linking to sp version: 1.2-5
# See what file types are supported by rgdal drivers
# ogrDrivers()$name
gdal.org
?readOGR
For more info on working with rgdal to load different types of spatial data in R see this excellent tutorial by Zev Ross.
Take a look at the file(s)
dir("data", pattern="sf_boundary")
## [1] "sf_boundary.dbf" "sf_boundary.prj" "sf_boundary.shp" "sf_boundary.shx"
sfboundary <- readOGR(dsn="data",layer="sf_boundary")
## OGR data source with driver: ESRI Shapefile
## Source: "data", layer: "sf_boundary"
## with 1 features
## It has 3 fields
# or
# sfboundary <- readOGR("data","sf_boundary")
# but not
#sfboundary <- readOGR(dsn="data/",layer="sf_boundary")
What type of sp object is sfboundary? How many features are in the object?
class(sfboundary)
str(sfboundary)
head(sfboundary@data)
Explore the object in the Envi window
sfboundaryHow?
sfboundaryplot(sfboundary)
sfboundary & sfhomes15_spplot(sfboundary)
points(sfhomes15_sp, col="red")
sfboundary & sfhomes15_spWhere are the points? What’s wrong?
plot(sfboundary)
points(sfhomes15_sp, col="red")
Compare the CRSs, are they the same?
proj4string(sfboundary)
proj4string(sfhomes15_sp)
proj4string(sfboundary)
## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
proj4string(sfhomes15_sp)
## [1] NA
sfboundary@bbox
## min max
## x 542696.6 556659.9
## y 4173563.7 4185088.6
sfhomes15_sp@bbox
## min max
## lon -122.51072 -122.36964
## lat 37.70854 37.80612
The #1 reason…
All sp objects should have a defined CRS
If not, one needs to be assigned to the object
This is called defining a projection
This doesn’t change the coordinates
All sp objects should have the same CRS.
This CRS should be defined in the proj4string slot.
When they don’t, they need to be transformed to a common CRS.
projecting or reprojection.Projection transformation returns a new spatial object with the transformed coordinates
GIS software - or software for working with geospatial data - will have a database of definitions for thousands of Earth referenced coordinate systems
and methods for using those CRS definitions to define or transform a CRS.
So, do we need to define the CRS of either
sfboundary or
sfhomes15_sp?
If yes, which one? How can you tell?
We need to know
sp objectbbox(sfhomes15_sp)
## min max
## lon -122.51072 -122.36964
## lat 37.70854 37.80612
Use the sp function CRS() to create a CRS object
Then, use the sp::proj4string() to assign it
# Define the CRS for sfhomes15_sp as WGS84
proj4string(sfhomes15_sp) <- CRS("+proj=longlat
+ellps=WGS84 +datum=WGS84 +no_defs")
The WGS84 proj4 string
CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
The WGS84 epsg code
CRS("+init=epsg:4326")
EPSG refers to the group that defined these
EPSG codes are more commonly used. Can you guess why?
# use an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326")
# or enter the parameter string
# proj4string(sfhomes15_sp) <- CRS("+proj=longlat
# +ellps=WGS84 +datum=WGS84 +no_defs")
Proj4 is the standard library for defining and transforming map projections and coordinate reference systems.
Here is an example file of proj4 CRS definitions
Once you set the CRS you can check it with the same function.
proj4string(sfhomes15_sp)
## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
proj4string() can get or set the CRS
CRS() is used to define the CRS of an sp object
The following are equivalent
# Define the CRS with the parameter string
proj4string(sfhomes15_sp) <- CRS("+proj=longlat
+ellps=WGS84 +datum=WGS84 +no_defs")
# Define the CRS with an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326")
epsg codes are used as short-hand for CRSs definitions
Are they both defined? Are they the same?
proj4string(sfboundary)
## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
proj4string(sfhomes15_sp)
## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
proj4string(sfboundary) == proj4string(sfhomes15_sp)
## [1] FALSE
Use sp function spTransform
Requires as input:
a sp object to transform with a defined CRS
a target CRS
Outputs a new spatial object with coordinate data in the target CRS
sfboundaryLet’s transform the CRS of sfboundary to be the same as that of sfhomes15_sp.
sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))
Why not the other way around?
Do they match now?
How will we know?
proj4string(sfhomes15_sp) == proj4string(sfboundary_lonlat)
## [1] FALSE
plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
Woo-hoo!
We can define the CRS of sp objects and transform them to the same CRS for mapping and spatial analysis.
Recap
# Define the CRS
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326")
#Transform the CRS
sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))
sfboundary_lonlatUse writeOGR to save sfboundary_lonlat to a new shapefile
See ?writeOGR for help
# write transformed data to a new shapefile
writeOGR(sfboundary_lonlat,
dsn = "data",
layer = "sfbounary_lonlat",
driver="ESRI Shapefile")
# is it there?
dir("data")
=
We want all data in the same CRS
Which one is best?
Geographic CRSs
4326 Geographic, WGS84 (default for lon/lat)
4269 Geographic, NAD83 (USA Fed agencies like Census)
Projected CRSs
5070 USA Contiguous Albers Equal Area Conic
3310 CA ALbers Equal Area
26910 UTM Zone 10, NAD83 (Northern Cal)
3857 Web Mercator (web maps)
See http://spatialreference.org/
Use this site to find EPSG codes and proj4 CRS strings
coordinates to convert the landmarks data frame to a SpatialPointsDataFrameproj4string to define it’s CRS as 3857 Web MercatorspTransform to transform the landmarks data to geographic coords (WGS84)# Convert the DF to a SPDF
coordinates(landmarks) <-c("X","Y")
# Define the CRS with an EPSG code for Web mercator
proj4string(landmarks) <- CRS("+init=epsg:3857")
# Transform the CRS to WGS84
landmarks_lonlat <- spTransform(landmarks, CRS("+init=epsg:4326"))
# map it
plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")
readOGR to load the shapefile “data/sf_highways.shp”class()proj4string()spTransform if need to project it to WGS84lineshighways <- readOGR(dsn="data", layer="sf_highways")
class(highways)
proj4string(highways)
highways_lonlat <- spTransform(highways, CRS("+init=epsg:4326"))
plot(sfboundary_lonlat)
lines(highways_lonlat, col="black")
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")
## OGR data source with driver: ESRI Shapefile
## Source: "data", layer: "sf_highways"
## with 246 features
## It has 5 fields
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
# QUESTIONS?
sp library - spatial objects and methods in RGDAL libary - readOGR, writeOGRSPDF